/*
 *  Solver.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 3/1/11.
 *  Copyright 2011 ICS - USI - Lugano. All rights reserved.
 *
 */

#ifndef NSE_SOLVER_Solver_H
#define NSE_SOLVER_Solver_H

#include "Array.h"
#include "Numerics.h"
#include "VectorHelper.h"
#include "CRS.h"

namespace nsesolver {

/**
 * \brief The Solver class ... TODO Doc
 * \tparam T 
 */ 
template< class Derived >
class Solver {
private:
	 inline const Derived& derived() const { return *static_cast< const Derived * >(this); }
public:
	
	template< class MatrixType >
	void step(const MatrixType &A, MatrixType &u, const MatrixType &f) const
	{
		derived().step(A, u, f);
	}
	
	
	template< class MatrixType >
	void stepSparse(const MatrixType &A, MatrixType &u, const MatrixType &f, const CRS< Array< int > > &crsIndex) const
	{
		derived().stepSparse(A, u, f, crsIndex);
	}
	
	
	template< class MatrixType >
	typename MatrixType::Scalar solveSparse(const MatrixType &A, MatrixType &u, const MatrixType &f) const
	{
		typedef typename MatrixType::Scalar Scalar;
		const Scalar err0 = 10000 + VectorHelper< Scalar >::norm(A*u - f);
		const Scalar tol = err0 * Numerics::eps< Scalar >();
		
		Scalar err = err0;
		int iteration = 0;
		
		CRS< Array< int > > crs(CRS< Array< int > >::CreateFrom(A));
		
		do {
			stepSparse(A, u, f, crs);
			if (hasToComputeError(iteration)) {
				Scalar prevErr = err;
				err = VectorHelper< Scalar >::norm(A*u - f);
				if (Numerics::equalsEps(err, prevErr)) {
					break;
				}
				
				std::cout <<"Iteration:" << iteration << ", err: " <<  err << ", tol" << tol << std::endl;
			}
			iteration++;
		}
		while (err > tol);
		
		return err;
	}
	
	
	template< class MatrixType >
	typename MatrixType::Scalar solve(const MatrixType &A, MatrixType &u, const MatrixType &f) const
	{
		typedef typename MatrixType::Scalar Scalar;
		const Scalar err0 = 10000 + VectorHelper< Scalar >::norm(A*u - f);
		const Scalar tol = err0 * Numerics::eps< Scalar >();

		Scalar err = err0;
		int iteration = 0;
		
		do {
			step(A, u, f);
			if (hasToComputeError(iteration)) {
				Scalar prevErr = err;
				err = VectorHelper< Scalar >::norm(A*u - f);
				if (Numerics::equalsEps(err, prevErr)) {
					break;
				}
				
				std::cout <<"Iteration:" << iteration << ", err: " <<  err << ", tol" << tol << std::endl;
			}
			iteration++;
		}
		while (err > tol);
		
		return err;
	}
	
	
	virtual
	bool hasToComputeError(const int iteration) const
	{
		return (iteration % 100) == 0;
	}
	
	
	
	/**
	 * The Destructor
	 */
	virtual
	~Solver()
	{
		//TODO
	}
	
};
}

#endif

